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ABSTRACT 

We study the persistence of warps in galactic discs in the presence of massive halos. 
A disc is approximated by a set of massive rings, while a halo is represented by a 
conventional n-body simulation. We confirm the conclusion of Nelson & Tremainc 
(1995) that a halo responds strongly to an embedded precessing disc. This response 
invalidates the approximations made in the derivation of classical 'modified tilt' modes. 
We show that the response of the halo causes the line of nodes of a disc that starts 
from a modified tilt mode to wind up within a few dynamical times. We explain 
this finding in terms of the probable spectrum of true normal modes of a combined 
disc-halo system. 
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1 INTRODUCTION 

At least 50 percent of spiral galaxies have warped discs, and 
it is likely that all galactic discs are slightly warped. Given 
the importance of the warp phenomenon and the fact that it 
is now nearly 40 years since the first warp, that of the Milky 
Way, was discovered, it is remarkable that there is still no 
generally accepted explanation of the phenomenon. Several 
recent reviews survey the literature of subject (e.g., Binney, 
1992; Nelson & Tremaine, 1996). 

To a first approximation most warped discs may be rep- 
resented by an ensemble of rigid, concentric rings. Each ring 
is spinning and coupled to the other rings and to non-disc 
material by gravity. The gravitational torques on the rings 
cause their spin angular momenta to precess, and the fun- 
damental problem posed by the phenomenon of warps is to 
understand how, despite this precession, the orientations of 
rings are coordinated, as they must be if the disc is to remain 
smooth and thin. 

An explanation that has enjoyed considerable popular- 
ity is that the rings are sufficiently tightly coupled to one 
another by gravity that they can precess together in the 
non-spherical gravitational potential of a massive halo as if 
they formed a rigid body (Toomre 1983; Dekel & Shlosman 
1983). Since in this model the height of the disc above the 
halo's equatorial plane varies sinusoidally in time, the disc is 
said to have an excited normal mode. The mode in question 
is a modification of a trivial zero-frequency solution of the 
linearized equations of motion of an isolated system of rings. 
This zero-frequency mode corresponds to a simple tilt of the 
system of rings with respect to the coordinate system em- 
ployed. Hence it is called the 'modified tilt' mode of the disc. 
The modified tilt modes of discs were thoroughly studied by 
Sparke & Casertano (1988). They found that some well- 
observed warps could be successfully modelled by modified 
tilt modes. 



Unfortunately, modified tilt modes are artificial in that 
they rest on the assumption that the halo potential, which is 
responsible for the precession of the disc, is unaffected by the 
disc's precession. In reality halo objects will tend to pick up 
energy from the time- varying potential of the disc, and any 
increase in the halo's energy is likely to be at the expense 
of the energy of the disc's modified tilt mode. Nelson & 
Tremaine (1995) estimated the rate of energy transfer from 
the disc to the halo and found it to be surprisingly large. 
They concluded that in realistic circumstances a modified 
tilt mode will be damped within one dynamical time of the 
disc. 

In this paper we use numerical simulations to test and 
clarify the semi-analytic work of Nelson & Tremaine. What 
they calculated was the rate at which halo stars would ac- 
quire energy if (a) the warped disc precessed precisely as pre- 
dicted by Sparke & Casertano, and (b) if the time-averaged 
potential governing the motion of halo particles is spheri- 
cal. The latter assumption is obviously unphysical and im- 
posed only in order to facilitate the calculations. However, 
intuitively one feels that the final result should not be sen- 
sitive to this assumption. By contrast the first assumption, 
that the disc precesses according to the prediction of Sparke 
& Casertano, is incompatible with the large energy trans- 
fer rate obtained by Nelson & Tremaine, since this demon- 
strates that effects neglected by Sparke & Casertano are in 
fact large. 

Simulations of disc-halo interaction are demanding nu- 
merically because the scales of the disc and halo are very 
different. Some numbers appropriate to the Milky way will 
illustrate this. The disc scale length is ~ 3kpc (Kent, Dame 
& Fazio, 1991) while the halo probably extends at least to 
the galactocentric distance r ~ 50kpc of the Magellanic 
Clouds (Fich & Tremaine, 1991). Moreover, it is essential 
for the halo mass and angular momentum to exceed those 
of the disc by a considerable factor. On the other hand, the 
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individual masses of halo particles should be smaller than 
those of disc particles in order that the halo's structure is 
sufficiently well defined in the neighbourhood of the disc. In 
particular, spurious stochastic acceleration of disc particles 
by halo particles must be avoided, as must significant two- 
body relaxation within the small fraction of the halo that lies 
within the disc. The difficulty of satisfying these conflicting 
requirements is illustrated by the simulations of Dubinski 
& Kuijken (f995). They simulated a disc-halo system with 
a direct n-body simulation in which 50 000 particles were 
assigned to the halo and 40 000 particles were assigned to 
the disc. Typically, their halo particles were more massive 
than their disc particles by a factor in excess of 6. In such 
a simulation the disc is inevitably so thick that one cannot 
meaningfully determine whether any warp's line of nodes is 
straight or has a tendency to wind up. 

Here we represent only the halo by an n-body simula- 
tion, while modelling the embedded disc by a series of rigid 
concentric rings. The latter approximation to disc dynamics 
was introduced by Hunter & Toomre (1969) and has since 
been widely used in theoretical work. However, it does not 
seem before to have been directly coupled to an n-body sim- 
ulation. It enables us to assign all our f 00 000 particles to 
the halo and to have a well defined line of nodes in a disc 
whose 100 individual rings have masses that are larger than 
those of halo particles by a factor of f 00. 

In Section 2 we set out the equation of motion of the 
disc and halo system. In Section 3, we describe numerical 
details and the initial conditions we have used for the halo 
and the disc. In Section 4 we check that our simulations 
can reproduce previous results and go on to show that a 
dynamical halo causes a warp to wind up rapidly when it 
starts from the configuration of the modified tilt mode for 
the case of a frozen halo. In Section 5 we explain this re- 
sult in terms of the probable normal-mode spectrum of the 
coupled disc-halo system. 



2 FORMALISM 
2.1 Disc dynamics 

The equation of motion of a disc star is 
A 2 z 9$ h 



(1) 



At 2 dz 

where "^h is the gravitational potential of the halo and / is 
the vertical force on the particle from the disc. For small 
excursions out of the xy-p\&ne, the halo term in (1) can be 
expanded as 
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F h = -d$ h /dz\ z= o 

k 2 . = d 2 $h/dz 2 \ z=0 - 

When we express the total time derivative in (I) in terms of 
partial derivatives we now have 
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where O(r) is the circular frequency at radius r. 

Following Hunter & Toomre (1969) we focus on the m = 
1 distortions of the disc by assuming that z is of the form 



z(t, r, 4>) = V2 [z c (t, r) cos(0) + z s (t, r) sin(0)] . 



(5) 



With this assumption the disc may be considered to consist 
of a nested sequence of rigid rings. Then \/2z c {t, r) is the 
height at which the ring of radius r passes over the posi- 
tive a>axis at time t, while s/2z s {t,r) is the height of the 
corresponding passage over the j/-axis. (The factors of \/2 
will ensure that the expression for the energy of the warp 
is natural.) When we similarly decompose Fh(t,r,(j>) and 
f(t,r,(f>) in the form (5), equation (4) yields coupled equa- 
tions of motion for z c and z B : 
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(6) 



We calculate the vertical forces between the rings to 
first order in z c and z s by splitting the force into two parts 

and (Toomre, private communication), f^ is the 
force on a ring that lies in the xy-plane from other, tilted 
rings, while is the force on a tilted ring when the other 
rings lie in the xy-pi&ne. For a ring of radius r, we have 
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where m r i is the mass of the ring of radius r' and 
A(r, r',9) = (r 2 + r' 2 + a 2 - 2rr cos 9) 3/2 . 



(7b) 



Here a is the softening of the force that arises because the 
disc has finite vertical thickness Az ~ a. 

If $d is the gravitational potential of the disc, the cir- 
cular frequency Q is given by 
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By analogy with the vertical forces, we decompose f2 d into 
contributions that are proportional to the two integrals that 
appear in equations (7a): 



n 2 d (r) = S (a) + S (b \ 
where 
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Hence we have finally that equations (6) may be written 
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Equations (11a) lead naturally to a definition of the 
warp's energy 



where 

if = (i„) 2 + {z si f . (13) 
and the two potential-energy terms are defined by 

N r 
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2.2 The halo 



We represent the halo by an ensemble of particles that move 
in the combined potential of the halo and disc. We expand 
this potential and the halo's density in spherical harmonics 
as follows, 



Ph(r, 9, <f>) = ^2 E p ™( cos ^ 

!=0 m=0 

x [Ai m {r) cos(m4>) + Bi m (r) sin(m0)], 
i 

1 = m=0 

x [Cim{r) cos(m0) + D im (r) sin(m^)]. 



(16) 



Here the pj™ are defined in terms of the conventional Legen- 
dre functions P ( m by 



_ (l-\m\)\ „ 

pi = VoTh)! P; 



(17) 



and the C'i m are related to Ai m by 
Ci m (r) = Cum{r) + C 2 i m (r) 
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An expression for the D tm in terms of the B tm is obtained 
by replacing C with P and A with B in these expressions. 

In the linear regime the vertical force on a ring is char- 
acterized by the force components, 
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and by the frequency, 
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where 
Pio = - 



r r 
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The halo's contribution to the circular frequency is 



«h = -^2Pi(0)E l0 (r). 



(22) 



2.3 The coupling between disc and halo 

The potential of a ring of radius r' may be written 

$(r,M)=E*^ P( ( COS7 )' 
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where the are coefficients to be determined, 
min(r', r) 



max r, r 



(23a) 



(23b) 



and 7 is the angle between the ring's symmetry axis and the 
direction (9,<f>). On the ring's symmetry axis, 7 = 0, we 
have 

t / „n Gm r i 
$(r,7 = 0) 



vV 2 + r 2 
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max r, r 
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Since P(l) = 1 for all I, the values of the follow immedi- 
ately from a comparison of equations (23a) and (24). 

The addition theorem for associated Legendre functions 
may be written 



P ; (cos 7) = ^2 -rPT '(cos 6 t )pT '(cos 6) 



(25a) 



x cos [m(0 — </>t)] , 



where 7 is the angle between the general directions (9t, 4>t) 
and (9, </>), and 



for m = 0, 
otherwise. 



(25b) 
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We set (8t, 4>t) equal to the polar angles of the ring's sym- 
metry axis 



e T (r ) = sin I -X I , 

4>t{t') = — + atan2 (z c , -z s ) , 

and eliminate P;(cos7) between equations (23a) and (25a 
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2.4 Energy integrals 

The potential energy of the system is V = Vhh + Vhd + Vdd, 
where, 



V hh = i / d 3 rp h $ 



\ hd = J d 3 rp h <i> d = V dh = J 



d rS d $ h , 



(28) 



Vdd= i / d 2 rS d $ d . 



The requirement that Vhd = Vdh provides a non-trivial check 
on the accuracy of the calculations. 

2.4-1 Self potential energy of the halo 

Vhh can be expressed in terms of the coefficients defined 
above of the spherical-harmonic expansions of p h and $ h : 

Vhh = ^2TkT 51 I ^™ AlmClm + eS rnBi m Di m ] , (29) 

1=0 m=0 

where e B m — 0, if m = and 1, otherwise. 

2.^.2 Interaction of the halo and the disc 

Similarly, the energy of interaction of disc and halo can be 
written 

2 1 f 

Vhd = ^ 2l - 1 ^ / drr 2 (e^A im Ci' m + e B m Bi m D' lm ) . (30) 

1=0 m = 

We decompose Vdh into two parts, a part V d u h associated 
with a flat disc that lies in the xy-p\&ne, and a part V d " 
associated with warping or tilting of the disc. We have 

Vdh= / drr d<^E d (r)<f> h (r,O,0). (31) 
J o Jo 



Setting 
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and using equation (16) we find 
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We calculate VJh by substituting equation (32) into the (28) 
and expanding $h in cylindrical polar coordinates (R, 4>, z): 
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On performing the integral over d 2 r we recover the first of 
equations (14). 

2.4-3 Self potential of the disc 

We similarly decompose the self potential energy of the into 
two parts, the energy V^ of an isolated flat disc, and the 
energy V d ™i associated with vertical displacements of the disc. 
Since V d u d is fixed, it may be neglected. V dd is given by the 
second of equations (14). 

2-4-4 Total energy 

The total energy is simply the sum of the warp energy, equa- 
tion (12), and the energy of the halo in the presence of an 
unperturbed disc: 

Etot = -Ewarp + -Ehalo, (35a) 

where 

£ halo = i ^ m ^ + Vhh + VS. (35b) 

i=l 

3 NUMERICAL DETAILS 

The coupled equations of motion of the rings and parti- 
cles were integrated with a symplectic leap-frog scheme that 
was developed according to the theory of Saha & Tremaine 
(1992). 

3.1 Representing the disc 

The disc is made up of 100 rings uniformly spaced in radius 
between n and rioo with rioo/fi = 18. Their masses mj are 
determined by the disc's surface-density profile £ d (r) an d 
the radii Ti of the rings. In most simulations we adopted 
the surface density profile of Sparke & Casertano (1988): 

M r<R t , 
S d (r) = E e- r/fld x i cos 2 (f ( ^ ~_^ t ) ) ) Rt < r < Ro (36) 

I 0, r> Rq. 

Thus interior to Rt the disc is exponential with scale length 
R d , while beyond R t its surface density tapers smoothly to 
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Figure 1. (a) The masses of the disc (solid curve) and the halo 
(dotted curve) inside a sphere of radius r. (b) The axis ratio c/a 
as a function of semi-major axis a. The open and full points are 
for the halo prior to and after insertion of the disc, respectively. 



The distribution of halo particles that has been pro- 
duced as described above is only in approximate dynamical 
equilibrium and has yet to be modified by the gravitational 
potential of the embedded disc. Consequently, we integrated 
the equations of motion of the halo particles for a time £ m 
during which the mass of a flat embedded disc is ramped 
from zero to its final value according to the formula 



M d (t) = \ [l - cos(7rt/i m )]M d (t m ), 



(39) 



where t m = 20/fl(R d ). Fig. la shows for our standard model 
the masses of the disc and the halo that lie inside a sphere of 
radius r. Fig. lb shows the axis ratio c/a of the isodensity 
surface in the halo that has semi-major axis a. 



3.3 Scalings 

Table 1 lists the numerical values of the parameters that appear 
in the defining equations of the disc and halo 

Our unit of time is determined by these choices: (i) n = 
r c = 1; (ii) the halo mass is unity and the disc mass is 0.1; (iii) 
Newton's constant G = 1. The circular frequency at R d is then 
Q(R d ) = 0.04 and the rms speed of halo particles is vo = 0.19. 
Consequently the halo's dynamical time is t dyn = ^Too/vq ~ 170. 



4 RESULTS OF THE SIMULATIONS 



3.2 The halo 

The bulge and halo are represented by N p = 10 5 particles. 
Since a general distribution function for a flattened halo will 
depend on a non-classical third integral for which no analytic 
expression is available, our halo model is constructed not 
from a distribution function but by allowing an approximate 
equilibrium configuration to relax towards an equilibrium. 
We construct an approximate equilibrium as follows. We 
start with a spherical system whose density profile varies 
approximately as r~ 2 . Specifically 



r/r t 



p{r) 



Po 




l+r 2 /r? 



for r < Too 
otherwise, 



(37) 



where r c is a suitable small radius, r t = 50r c and = 65r c . 
The isotropic distribution function f(E) of this system is 
obtained from Eddington's formula [eq. (4-140a) of Binney 
& Tremaine, 1987] and sampled by standard procedures. 

This spherical system is next flattened by multiplying 
each z coordinate by q < 1. We now multiply the z- velocity 
of each particle by a factor a and multiply all x- and y- 
velocities by a factor j3 to ensure that the tensor virial theo- 
rem is satisfied. Specifically we require that the components 
Tu of the kinetic-energy satisfy 2Tu — —Wu, where Wa 
is a component of Chandrasekhar's potential-energy tensor. 
Since T zz oc o? and T xx = T yv oc /3 2 , one easily finds that 



W zz 

W XX + Wyy 

4T° ' 



(38) 



where T° z is any non-vanishing component of the kinetic- 
energy tensor of the original spherical system. 



4.1 The case of a static halo 

We first check that our rings are capable of reproducing the results 
of Sparke & Casertano (1988) by freezing the configuration of the 
halo once the disc mass has been ramped up to its full value. 
Equation (5) of Sparke & Casertano was solved for the pattern 
speed Q p and vertical displacement z(r) of the modified tilt mode 
of our disc in this fixed halo. The equations of motion (11a) of 
the disc were then integrated from this initial condition. Fig. 2 
shows that the numerical solution is exactly as it should be for a 
normal mode: the disc maintains it shape to high accuracy while 
precessing at the derived pattern speed. Fig. 3 shows that the 
warp energy is conserved to high accuracy - it decreases by 5 
parts in 10 5 during 13 revolutions of the ring at r = R d . 

These results demonstrate that we have accurately modelled 
the disc-halo interaction and are correctly following the dynamics 
of the disc. 



4.2 Forced precession of the disc 

Next we check the calculation of Nelson & Tremaine (1995) by 
evolving the halo while the disc is caused to precess at a steady 
rate in a given warped configuration. The latter was the ap- 
proximate normal-mode configuration considered by Nelson & 
Tremaine. 

Fig. 4 shows the consequent evolution of the halo energy that 
is defined by equation (36). The disc's warp energy is -E W arp — 
3 X 10~ 4 , so, in so far as the straight line of slope = 2.8 X 10 — 7 
that is plotted in Fig. 5 provides an adequate approximation to 
the dependence of -Ehalo on time, the decay rate of the warp is 
T = -B halo /£ warp = 9 X 10~ 4 and T/n(R d ) = 0.02. This value 
of T lies exactly in the middle of the range of values for the same 
quantity that Nelson & Tremaine plot in their Fig. 2. 

This calculation suggests both that we are integrating the 
equations of motion of the halo particles accurately, and that the 
approximations employed by Nelson and Tremaine to estimate 
the effect of the warp on the halo are not misleading. 
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Tabic 1. Parameters 
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Figure 2. Plots of z(t, r,cf> = Q p t) for t = (curve) and t = 2000 
(triangles) when the disc is allowed to evolve from the configura- 
tion of a normal-mode for time t in a frozen halo (2000f2 p = 164°). 
The open squares show z(t, r,<f> = S7 p t + ^tt) at t = 2000. (The 
corresponding values of z for t = arc identically zero.) 



4.3 A live disk in a live halo 

Fig. 5 shows the evolution of the curve of nodes of live disc that 
evolves in a live halo from an initial condition that is a normal- 
mode in the sense of Sparke & Casertano. It is apparent that the 
warp winds up within a few dynamical times. Fig. 6 shows, as a 
function of time, the rms vertical displacement \z\ = (zj? + z 2 ) 1 ^ 2 
for this simulation. This does not evolve dramatically. Only in 
the innermost third of disc is there any tendency for the rms 
displacement to decrease with time. In the outermost quarter of 
the disk the rms displacement clearly increases slightly through 
the duration of the simulation. It turns out that the growth in 
\z\ at large radii predominates in the sense that _E W arp tends to 
increase. 
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Figure 3. Warp energy [equation (12)] as a function of time for 
a disc evolving in a frozen halo. 
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5 CONCLUSIONS 

We have used a program that uses rings to represent the disc and 
particles to represent the halo to model the dynamics of a disc 
that is embedded in a live massive halo. We have checked the 
correctness of the program as regards the disc and its interaction 
with the halo by showing that when the halo is frozen, and the 
disc is started in a Sparke— Casertano normal mode, it precesses 
rigidly as the normal-mode calculations require. We have checked 
the correctness of the program as regards the halo by showing that 
when the disc is rigid and forced to precess at a constant rate, the 
halo acquires energy at the rate predicted by Nelson & Tremainc 
(1995). 

When the halo is live and the disc is started from the con- 
figuration of a normal mode, the warp is not rapidly damped as 



Figure 4. The energy of the halo [equation (36)] when a warped 
disk is forced to precess at a constant rate. The dotted line is 
the fit used in the text. Notice that Swarp/^halo ~ 0.02, so 
calculating energies to the required accuracy is computationally 
challenging. 

Nelson & Tremaine predicted, but winds up within a few dynam- 
ical times, whilst retaining or even enhancing its 'warp energy'. 
How may this result be understood? 

We have to consider two distinct dynamical systems: (a) 
the disc in a frozen halo, and (b) the disc in a live halo. Let us 
call (a) the Frozen System and (b) the Live System. The tight 
coupling between the disc and halo that is demonstrated by the 
calculation of Nelson & Tremaine already implies that the normal 
modes of the live system are not close to the normal modes of the 
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Figure 5. A warp in a live halo rapidly winds up even though its initial configuration would be a normal mode if the halo were static. 
Here we plot at several times the curve of nodes that is defined by <f>(r) = n + atan2[z c (r), — z B (r)] for a live disc that evolved in a live 
halo from a configuration that was a normal mode in the sense of Sparkc & Cascrtano. 



frozen system. Hence, if we were to express the configuration that 
corresponds to any low-order normal mode of the frozen system 
as a linear combination of normal modes of the live system, the 
sum would not be dominated by a single term. The configuration 
described by the sum would evolve significantly on a timcscalc 
that is inversely proportional to the biggest frequency difference 
Au) max between the normal modes that have non-negligible am- 
plitudes in the sum. The rapidity with which a warp of the live 
system winds up from a normal-mode configuration of the frozen 
system implies that Aaj max is large. 

To understand why this should be so, consider the probable 
form of the normal modes of the live system. At any radius the 
density of the halo has a strong tendency to peak at the vertical 
location of the disc, because the latter is extremely dense. Equiv- 
alcntly, a substantial input of energy is required to displace the 
disc from the halo's local equator. By contrast, relatively little 
energy is required to tip the outer disc and halo with respect to 
the inner disc and halo - the outer and inner parts of a flattened 
galaxy are kept in angular alignment principally by the interac- 
tion between the mass quadrupole of the outer galaxy with the 
inner galaxy's potential quadrupole, which is attenuated by the 
cube of the ratio of the characteristic radii of the inner and outer 
parts of the galaxy. 

Now it is generally true that deformations of a system that 
are associated with small energy increments are associated with 
low-frequency normal modes, while those that are associated with 



large energy increments are associated with high-frequency nor- 
mal modes. So a normal mode that displaces the disk from the 
local equator of the halo will have a much higher frequency than 
one that tips the outer parts of the galaxy with respect to the 
inner parts. 

To generate a normal mode of the frozen system we have 
to do two things: First, we generate the warp in the disc by 
rotating the outer galaxy with respect to the inner galaxy. Then, 
we eliminate the twist just introduced into the halo by displacing 
the halo with respect to the disc at both large and small radii. 
Thus, when expressed as a sum over normal modes of the live 
system, a normal mode of the frozen system involves both high- 
frequency and low-frequency terms in an essential way. From this 
it follows that it will wind up rapidly, essentially at the frequency 
of the high-frequency normal modes. 

The argument we have just given explains why our initial 
warped configurations rapidly wind up. It also suggests that true 
normal modes of the full live system may involve configurations 
of the disc that are very similar to those of Sparke-Casertano 
normal modes. It is in their halo configurations that these true 
normal modes would differ essentially from Sparke-Casertano nor- 
mal modes. It is likely that these differences in halo structure will 
give rise to the true normal modes having significantly lower fre- 
quencies than do Sparke-Casertano normal modes. We hope soon 
to return to this prediction. 
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Figure 6. The rms vertical displacement |z| = (zj? + z 2 ) 1 / 2 as a function of radius for the simulation whose curves of nodes are shown 
in Fig. 5. 
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